Variational approach for Stokes flow through a two-dimensional non-uniform channel

A variational approach is proposed to study the Stokes flow in a two-dimensional non-uniform channel. By using the stationarity of the Lagrangian, the Euler-Lagrange equations are established which leads to a simple set of ordinary differential equations to provide an estimate for the average pressure drop explicitly in terms of the channel shape function. The results for the pressure drop show an excellent agreement with the second-order extended lubrication theory. A higher-order formulation further improves the accuracy of the results for the pressure drop along the channel.

computational data, which further highlights the upgraded model's interest in precisely capturing the velocity profiles and pressure drop along a channel.Luca & Llewellyn Smith 31 solved Stokes-flow equations in a two dimensional channel with a linear expansion using the Unified Transformation Method (UTM) by splitting the domain into three convex sub-domains, providing quasi-analytical solutions.The results obtained using the UTM for the pressure drop were compared to those obtained using the ELT adapting the analysis of Tavakol et al. 30 .Hinojosa et al. 32 recently presented an extended lubrication approximation for arbitrary shapes while the work of Tavakol et al. 30 was restricted to asymmetric channels.Few more studies on lubrication approximation for Newtonian and viscoelastic rheology were reported recently by Christov et al. 33 , Dewangan & Datta 34 , Boyko & Stone 35 and in the references there.
Variational principles often play an important role in basic physics to provide a deeper understanding of natural laws.In research based on variational principles, there are two groups of approaches: the first one finds the minimal entropy production for dissipative systems, while the second one is created from the minimum action integral in classical mechanics of energy-conserving systems.Fermat's principle, which is analogous to Snell's rule of light refraction and was put forth in the 17th century to anticipate the path of light, is where the minimum action integral got its start.Hamilton and other scientists later extended Maupertuis' concept of the minimum action integral in mechanics, which was based on this theory.The second set of variational principles originates from the work of Prigogine and his associates 36 , wherein systems with thermal nonequilibrium are crucially dependent on the condition of minimum entropy production.A variational principle is applied to the problem of instability in the context of thermal convection [37][38][39] .The equilibrium between the work performed by the buoyancy force and the viscous dissipation in the convective flow is the critical condition of the Benard convection, and the critical Rayleigh number is determined by the stationary value of the ratio of these quantities 37 .Takaki 40 used the variational principle, which necessitates the stationarity of the difference between the work performed by the pressure force and the entropy production, to construct the Stokes equations and the continuity equation.Wang et al. 41 established the form of the Lagrangian function using the variational principle.Other recent use of variational approach can be found in the works of Xu et al. 42 and Zhou & Doi 43 .
In this paper, we focus on the Stokes flow in a two-dimensional non-uniform channel, as shown schematically in Fig. 1.The problem is analyzed by constructing a functional, named Lagrangian for Stokes flow using the variational principle proposed by Takaki 40 and was implemented by Wang et al. 41 .Considering the stationarity of the Lagrangian with respect to the set of unknown variables, the Euler-Lagrange equation which estimates the average pressure drop between the two ends of the non-unform channel is formulated.We compare our results with those obtained from the second-order ELT presented by Tavakol et al. 30 .We also perform finite volume numerical simulations based on the solution of Stokes equations to validate our results.
This paper is organized as follows.In the section 'Formulation of the problem and mathematical model' , the problem of Stokes flow in a two-dimensional non-uniform channel is formulated using the variational approach.Section 'Numerical simulations' describes the details of the numerical method used here.The results obtained in terms of variation of the pressure drop for a constricting channel are presented in the section 'Variation of pressure drop for a constricting channel' .The results for the pressure drop variation in an expanding channel are described in section 'Variation of the pressure drop for an expanding channel' .Section 'Estimate of the improvement of the solutions' presents the estimation and prediction of the improved results obtained in terms of the pressure drop using the velocity profile proposed by Hinojosa et al. 32 .The salient conclusions of the paper are summarised in section 'Conclusions' .

Formulation of the problem and mathematical model
We consider a steady, two-dimensional pressure-driven flow of an incompressible Newtonian fluid in a channel of length 2L 0 bounded by a planar horizontal wall on one side and a corrugated wall on the other.The material properties of the liquid are density ρ , dynamic viscosity µ and the kinematic viscosity ν = µ/ρ .The x-axis of the chosen reference frame coincides with the planar wall and the y-axis is normal to the x-axis and points into the www.nature.com/scientificreports/fluid.The location of the channel walls are given by y = 0 and y = h(x) = h 0 H(X) , where h 0 is a characteristic channel height and H(X) is the dimensionless channel height.The wall shape function h(x) is assumed to be at least C 2 continuous in the entire stretch of the channel.We also assume that the Reynolds number R = Q/ν based on the volumetric flow rate through the channel Q is small, so the flow is governed by Stokes equations where u = (u, v) is the velocity field, p is the pressure and µ is the fluid dynamic viscosity.
The Classical Lubrication Approximation (CLA) yields 30,32 which can be simplified as where û(x) = Q/h(x) and g(x) are functions of x which are to be determined.Following Takaki 40 and Wang et al. 41 we construct the Lagrangian for Stokes flow in the corrugated domain as Using the expressions (5)- (7) for the CLA of Stokes flow, we obtain Combining the above integrals we obtain the Lagrangian L as a functional of û and g in the form Using the methods of calculus of variations we obtain a set of two Euler-Lagrange equations as (1) Vol:.( 1234567890) www.nature.com/scientificreports/which reduce to respectively.The first Euler-Lagrange Eq. ( 15) results in which is automatically satisfied since û(x)h(x) = Q is a constant flow rate Q along the channel, and therefore Substituting Eqs. ( 17) and ( 18) into the second Euler-Lagrange Eq. ( 16) we obtain Hence, it is possible to obtain an expression for the average-pressure gradient p′ (x) along the channel as In order to obtain the dimensionless average pressure gradient, the following dimensionless variables are introduced In terms of these non-dimensional variables, we obtain the dimensionless pressure gradient as where is the aspect ratio of the channel.Thus, the pressure drop P can be determined as

Numerical simulations
We compare our results for the pressure drop obtained analytically with the corresponding results obtained using direct numerical simulations.To do this, we solve numerically the full time-dependent Navier-Stokes equations for a two-dimensional incompressible flow in a channel with prescribed boundary conditions.As a result of this solution starting with an arbitrary chosen initial condition, we look for a steady solution of the system.The continuity and Navier-Stokes equations (in the Stokes theory limit for the exceedingly small Reynolds number) are brought into dimensionless form using the scaling and dimensionless parameters introduced in Eqs. ( 22) and ( 24) to obtain www.nature.com/scientificreports/ In order to convert the corrugated geometry of the channel shape into a conventional rectangular one a mapping 44,45 is applied to Eqs. ( 26)- (28).Under the mapping (29), Eqs.(26-28) are recast as In our direct numerical simulations, a staggered-grid-based finite volume method (FVM) is employed to solve the set of Eqs. ( 30)- (32).By progressing the flow field variables via a series of lesser time steps of size 0.001, a time-marching numerical approach is obtained.It is found that the flow field achieves a steady state independent of the stipulated initial conditions after a brief transient for the range of parameter values studied here.A cyclic pressure correction approach, the SIMPLE algorithm, is used to solve the discretized Eqs. ( 30)-( 32) [44][45][46] .
The discretized continuity equation is transformed into a Poisson equation for the pressure correction, which is then solved iteratively using the SOR (Successive over relaxation) scheme until the desired accuracy is achieved.This creates the pressure link between the momentum and continuity equations.The convergence criterion is defined as i,j |< � , where � = (U, V , P) , = 10 −6 , n denotes the iteration number, and i, j stand for the number of the node of the computational grid.The numerical simulations are performed in the non-dimensional transformed domain (η ∈ [−1, 1]) × (ξ ∈ [0, 1]) with the non-uniform spacing with denser grid points in the wall vicinity to capture the gradients of the velocity components.The typical grid spacing used in our simulations in the η direction is varied between 0.005 and 0.01, whereas the grid spacing in the ξ direc- tion is uniform in the range between 0.01 and 0.02.The dimensionless channel height is considered H 0 = 1 .Additionally to a no-slip, no-penetration conditions on the walls, the boundary conditions are defined using the Poiseuille velocity profile at the inlet, and a vanishing outlet gauge pressure 30 .Thus, the boundary conditions in the transformed region [−1, 1] × [0, 1] takes the form

Results and discussion
In this section, we present our analytical results for the pressure drop based on the variational formulation derived in Section Formulation of the problem and mathematical model, for a non-uniform channel.We also aim to compare our results with those of the CLT (Classical Lubrication Theory) and a higher-order ELT (Extended Lubrication Theory) developed by Tavakol et al. 30 .We also validate the predictions of our variational model by performing 2D numerical simulations using the finite-volume approach described in Section Numerical simulations.
To validate and compare our results derived by the variational approach, we choose the same example considered by Tavakol et al. 30 where the normalized wall shape function is given by where represents the dimensionless amplitude of the channel height variation; > 0 and < 0 stand for the cases of constricting and expanding channel shapes, respectively.For the purpose of our fully numerical solution only we extend the channel by www.nature.com/scientificreports/For a fixed value of the flow rate Q, the dimensionless pressure drop between X = −1 and X = 1 in the case of the shape function (34) is obtained by Tavakol et al. 30 where P 0 , P 2 and P 4 represent the pressure drop at leading order in δ , and the second-and fourth-order lubrication approximation, respectively.The expressions for P 0 , P 2 and P 4 are explicitly given in 30

Variation of pressure drop for a constricting channel ( > 0)
In this section, we explore the effect of the dimensionless amplitude and the aspect ratio δ on the pressure drop in the case of the channel with a concave shape.We compare the values of the pressure drop derived using the variational approach with those obtained for the extended lubrication theory (ELT) of various orders 30 and our simulation results.
First, in Fig. 2 we present the non-dimensional pressure drop �P = �p/(µQL 0 /h 3 0 ) between the two ends of the channel as a function of with a fixed δ = 1 .Our variational results along with those for various-order ELT together with our simulation results are all presented in Fig. 2. It is evident from this Fig. 2 that the pressure drop increases rapidly with an increase in , as expected, since the flow rate is constant through any cross section along the channel.The most important observation from Fig. 2 is that our results for the pressure drop derived from the variational approach agree very well with the second-order ELT for a wide range of with δ = 1 .However, it is also evident from Fig. 2 that the experimental results 30 are in a better agreement with the second-order ELT (alternatively, with our variational results) rather than with those of the fourth-order ELT.A sample of our simulation results for the velocity magnitude with = 0.5 and δ = 1 are presented in Fig. 3.
In Figs.4a and b, we compare between our analytical results, the results of various-order ELT solutions and those of our numerical simulations.It is evident from Fig. 4a that when δ is small, e.g.δ = 0.5 , the difference between the second-and fourth-order ELT results is negligible and they fit our results arising from our both variational and numerical simulation results very accurately.For higher values of δ , say for δ = 1.5 , one can observe that the results derived from the variational approach agree with those of the second-order ELT.It is important to mention that the non-dimensional pressure drop derived from the variational approach and various-order ELT deviate from the simulation results above a certain value of (say, = 0.5).Comparison between the results for the dimensionless pressure drop obtained for δ = 1 using our variational approach based on the parabolic velocity profile (the dotted curve), our numerical finite-volume simulations (squares ), and the second-( + ) and fourth-order ELT (the solid curve) derived by Tavakol et al. 30 for a constricting channel ( > 0 ).The triangular symbols represent the experimental results obtained by Tavakol et al. 30 .

Variation of the pressure drop for an expanding channel ( < 0)
This subsection is devoted to discussion of our results based on the variational formulation and to their comparison with various-order ELTs 30 and with our numerical results for an expanding channel of a convex shape.Such comparison is depicted in Fig. 5 for a channel defined with a convex shape function ( < 0 ) for δ = 1 .The convexity of the shape alters the characteristics of the p − dependence, as shown in Fig. 5.As increases, the pressure drop between the two ends reduces to conserve the flow rate in an arbitrary cross section of the channel.It is evident from Fig. 5 that similar to the results obtained for contracting geometries discussed above, our variational results for the dimensionless pressure drop agree very well with the second-order ELT in case of an expanding channel for a wide range of with δ = 1 .However, for a channel with a sharp convexity, the results obtained from the fourth-order ELT show a better agreement with our numerical simulations.Figure 6 depicts the non-dimensional pressure drop p between the two ends of an expanding channel as a function of for δ = 0.5 and δ = 1.5 in its left and right panels, respectively.It is observed in the left panel of Fig. 6 that for a small value of δ , e.g., δ = 0.5 , the differences between the results of our numerical simulations, and those of the second-and fourth-order ELT is very small and they agree with our variational results very well (within 4% error).For a higher value of δ , e.g., δ = 1.5 , the results are presented in the right panel of Fig. 6.The comparison reveals that our variational results for the pressure drop agree very well with the second-order ELT for a broad range of , even for δ > 1 .However, our simulation results agree with our variational results (alternatively, with the second-order ELT) up to a small value of , 0 > > −0.3 (within 5% accuracy), and the discrepancy begins to increase rapidly with an increase in the channel convexity.It is also observed from the right panel of Fig. 6 that in the case of an expanding channel with δ = 1.5 , the fourth-order ELT can not provide a better estimate than the second-order ELT in the range of moderate .
Finally, one should bear in mind that the case of an expanding channel represents, in general, a more complex setting than that of a constricting one.The possibility of a flow separation may invalidate the variational approach based on a parabolic velocity profile.Figure 7 shows the evidence of flow separation near the wall crest of the channel for an expanding channel with = −0.60 and δ = 1, as obtained by our fully numerical solution of the flow.Flow reversal magnified in the inset to Fig. 7 for a better view causes a significant difference between the fully numerical results and our results based on the variational approach for a higher value of | | .It is worth mentioning that for a low aspect ratio of the channel, e.g., δ = 0.5 , flow separation does not take place, and thus our variational results are expected to show a good agreement with fully numerical results, as seen in Fig. 6a.Variation of the dimensionless pressure drop with as obtained using the variational approach based on the parabolic profile (the dashed curve), our finite-volume simulations ( ), and second-( + ) and fourth- order (the solid curve) lubrication approximations both derived by Tavakol et al. 30 for an expanding channel ( < 0 ).Left panel: δ = 0.5 ; Right panel: δ = 1.5.Comparison of the results obtained for δ = 1 for the dimensionless pressure drop using our variational approach based on the parabolic velocity profile (the dotted curve), our numerical finite-volume simulation (squares ), and the second -( + ) and fourth-order (the dashed line) ELT derived by Tavakol et al. 30 .
The solid curve represents the results obtained using our variational approach based on the velocity profile derived by Hinojosa et al. 32 for a constricting channel.The inset allows a better view of the five different sets of results in a smaller domain of .

Figure 1 .
Figure 1.Sketch of the problem.

Figure 7 .Figure 8 .
Figure 7. Streamlines map obtained using finite-volume simulations for = −0.60 and δ = 1 .The inset shows a magnifying view of the streamlines near the reverse flow region. 2dX.